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Abstract. It is known that symmetries inherent in the statistical theory of ideal (i.e., non- 
dissipative) magnetohydrodynamic (MHD) turbulence are broken dynamically to produce coherent 
structure. Previous numerical investigations are extended to study decaying MHD turbulence. Here, 
we find that coherent structure also arises and persists in the presence of dissipation. Thus, the mag- 
netic dynamo may be due to a dynamically broken symmetry that is inherent to MHD turbulence. 
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INTRODUCTION 

In this work, incompressible, homogeneous 3-D MHD turbulence is examined in terms 
of truncated Fourier models. Physical fields such as the turbulent magnetic induction b 
are represented by truncated Fourier series on an TV 3 grid: 

b(x) = £ b(k) e ikx . (1) 

0<|k|</f 

The wave vectors k in k-space have integer components between —TV/2 and +TV/2, 
while the corresponding discrete set of position vectors x in x-space have components 
that are integer multiples of 2n/N, where TV is the number of grid points in each 
dimension. The ‘isotropic truncation’ radius is K, which is chosen to ensure that retained 
coefficients fit within a sphere in k-space; the exact value of K is dictated by de-aliasing 
requirements [Patterson 1971]. 

These models place MHD turbulence in a periodic box that can represent both a small 
volume interior to a larger volume containing plasma, and also a complete mass of 
plasma within a 3-torus (such as the universe). Furthermore, isolated plasma (such as 
a star) may be thought of as residing in box large enough so that any plasma flows 
or magnetic fields outside of or on the surface of the box have negligible energy. For 
this later case, a Fourier model can serve as an approximate surrogate. Thus, Fourier 
models are of value by providing at least qualitative information about geophysical and 
astrophysical turbulence, especially if dissipation is included. 

Statistical studies of such Fourier models were initiated by T. D. Lee, in a prescient 
paper [Lee 1952] that demonstrated the existence of a canonical ensembles based on 
the ideal invariants in the phase space defined by independent Fourier modes. The 
total energy E is one ideal invariant, and two other invariants were soon discovered: 
magnetic helicity H M [Elsasser 1956] and cross helicity He [Woltjer 1958]. Once com- 
putational resources were sufficient, the basic ‘absolute equilibrium ensemble theory’ 



was described for three-dimensional (3-D) MHD turbulence [Frisch 1975]. The bro- 
ken symmetry and non-ergodicity contained in this theory was subsequently recognized 
[Shebalin 1982, 1983, 1989, 1994, 1996, 1998, 2002], 

A recent paper [Shebalin 2006] examines ideal, homogeneous, incompressible, 3-D 
MHD turbulence in a rotating frame of reference with and without a mean magnetic 
field. There turn out to be five general cases of incompressible, homogeneous 3-D 
MHD turbulence, differentiated by the absence or presence of a mean magnetic field 
B 0 and/or an overall system rotation Q 0 ; these cases are shown in Table 1. There, a 
new invariant, the parallel helicity Hp = H (: — oH M is indicated, which appears when 
Q 0 = ctB 0 ^ 0. (For further details, including details of very long-time 32 3 simulations, 
please see [Shebalin 2006].) 

Here, dissipation is introduced through non-zero viscosity v in eq. (2) and magnetic 
diffusivity rj in eq. (3). Results from a long-time 64 3 simulation are presented and large- 
scale coherent structures, as seen in previous ideal runs, are still present. In effect, this is 
a ‘magnetic dynamo process’ that is inherent to MHD turbulence, a process that requires 
an initial energetic stirring of a magneto-fluid, but does not require any continuous 
external forcing, such as convection. A magnetic dynamo may thus be related more 
to an inherent, dynamically broken symmetry, than to external forcing. 


BASIC THEORY 

The well-known non-dimensional form of the MHD equations in a rotating frame of 
reference with constant angular velocity Q 0 and a mean (i.e., uniform and constant) 
magnetic induction B 0 are (where V • u = V • b = 0, (D = Vxu and j = V x b) 

= V x [u x (co + 2Q 0 ) +j x (b + Bo)] + vV 2 co, (2) 

ot 

Tu 

— = Vx[ux(b + B 0 )] + t]V 2 b. (3) 

dt 

In equations (2) and (3), u, ft), b and j are the turbulent velocity, vorticity, magnetic 
induction and electric current density, respectively, of the magneto-fluid. 

If we set v = 0 in (2) and rj = 0 in (3), we obtain the equations of ideal MHD 
turbulence. While any real flows have v ^0 and )] / 0, a study of ideal turbulence 
produces some interesting theoretical results that may help in understanding aspects of 
real turbulence, particularly at larger scales of the flow where dissipation is minimal. 


IDEAL MHD TURBULENCE 

Define the volume average of a quantity (f> (r) in a periodic box of side length 2k as 

[<p] = (2 n)- 3 j (j)(r)d 3 x. 


(4) 



TABLE 1. Invariants for ideal MHD turbulence. 


Case 

Mean Field 

Rotation 

Invariants 

I 

B 0 = 0 

£2 0 = 0 

E, H c , H m 

II 

B o ^0 

£2 0 = 0 

E,H C 

III 

B 0 = 0 


e,h m 

IV 

B o ^0 

Q . 0 x B 0 y 0 

E 

V 

B 0 yo 

o 

II 

o 

a 

E,H P 


Since all functions are periodic, we can use equations (2) through (4), along with 
integration by parts to derive the following relations (see [Shebalin 2002]): 


dE 

dt 

dH c 

dt 

dHf/i 

dt 


= —2 (vQ + r/J), 


= £2 0 • [b x u] — - (v + q) [j • o] , 


= B 0 ■ [b x u] — rj [j • b . 


(5) 

(6) 

(7) 


Above, we have the (volume-averaged) energy E, enstrophy Q, mean-squared current J. 
cross helicity Hk and magnetic helicity Hm'- 

E='-[u 2 + b 2 ], Q=i[|»| 2 ]. r=i[|j| 2 ], ffc=l[u-b], Hm — ( [a-b] . (8) 


At this point we note that if Q 0 = <jB 0 , i. e . , if Q 0 and B 0 are non-zero and parallel, 
then eq. (6) can be added to —o times eq. (7) to yield 


^ ( v + *7 ) U ' ®] + a7 7 U ' b] , H P = H C - oH m . ( 9 ) 

The quantity ELp will be called the ‘parallel helicity.’ 

When v = rj = 0, equations (5), (6), (7) and (9) lead us immediately to the invariant 
integrals for MHD turbulence for various values of B 0 and Q 0 . These various cases are 
placed in Table 1 for easier reference. We see that as a mean field is applied or an overall 
rotation is imposed or both, the number of invariants drop from three to two to one. 

The probability distribution function D of the absolute equilibrium ensemble, for Case 
I in Table 1, has the general form 


D ~ exp(—aE — f5Hc — yH M ). 


( 10 ) 


The function D is normalized by requiring that its integral over the whole phase space 
defined by the independent Fourier coefficients be unity. In regard to the other cases in 
Table 1, eq. (10) applies as follows: Case II, set y = 0; Case III, set /3 = 0; Case IV, set 
/3 = y=0; and in Case V, set y= — a/3, so that eq. (10) becomes D ~ exp(—aE — f5Hp). 

Thus, to determine various ensemble predictions for Cases II through V in Table 1, 
we merely take previously determined results for Case I [Frisch 1975], [Shebalin 1989, 
1994, 2002, 2006], [Stribling 1990] and substitute for /3 and y as indicated above. 



(a) (b) 



FIGURE 1. Phase portraits for Case I in Table 1. Re vs Im parts of b(k), where k ■ b(k) = 0, are plotted 
for (a) k = (1,0,0) and (b) k = (0,1,1), from t = 0 to 300. The runs were on a 64 3 grid with v = t] = 1/250. 
The initial points ( t = 0) are indicated by o. Run 1 A results are in gray and those of Run 1 are in black. 


Using (10), ensemble predictions for moments of the u(k) and b(k) can be made for 
the different cases in Table 1. Denoting the ensemble prediction of a quantity <2(k) by 
(<2(k)), the predicted first-order moments are (u(k)) = (b(k)) = 0 for all cases. The 
ensemble predictions for second-order moments are not needed here, but are given in 
[Shebalin 2006] for the five cases listed in Table 1. 

An apparent challenge in using eq. (10) is that Hq , H M and H P are pseudoscalars 
under the discrete symmetries of charge inversion C and spatial inversion (or parity) P, 
in which case D in eq. (10) may not be invariant under C and P. However, it can be shown 
that the inverse temperatures /3 and y in eq. (10) are also pseudoscalars under C and P 
[Shebalin 1996, 1998, 2002, 2006], Thus, the probability density function indicated by 
eq. (10) is symmetric under C and P , as are the MHD equations (2) and (3). 

The basic result is that when eq. (10) is used for ensemble averaging, the integral 
over phase space includes disjoint components identified by plus and minus signs of the 
helicity. The number of disjoint components for the cases in Table 1 is 2 M , where M 
is the number of invariant helicities listed in each of the 5 cases (i.e., Case I has four 
components, Case IV has one component, and Cases II, III and IV have two components 
each). This disjoint structure was first pointed out in [Shebalin 1998]; applying the 
Birkhoff-Khinchin theorem [Khinchin 1949] tells us that only Case IV is ergodic, while 
Cases I, II, III and V are non-ergodic. 

However, when a numerical simulation is run or a physical system evolves, we can 
start with either the plus or the minus sign of a given helicity, but not both. Thus the 
symmetry of the theory is dynamically broken, and mean values of the Fourier coeffi- 
cients with respect to time are not necessarily zero, even though this is the ensemble 
prediction. Numerical examples of this are presented next. 
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FIGURE 2. Energy in modes with given k 2 < 4 vs time. E tot is the energy in all modes. 


REAL & IDEAL NUMERICAL RESULTS 

A Fourier spectral transform method [Patterson 1971] with 3 rd -order time-integration 
[Gadzag 1976] was used for numerical studies. Ergodicity can be examined using phase 
diagrams, i. e. , plots of the real versus the imaginary parts of a Fourier coefficient. 
Consider the time history of several Fourier modes b(k) from two 64 3 Case I runs 
(B 0 = Q 0 = 0, t - 0 to 300): Run 1 A (v = T] = 0) and Run 1 (v = 7] = 0.004); both runs 
had identical initial conditions: E(k ) ~ k 4 exp(0.5k/k p ), k p m 6, with random phase. 
Initially, total kinetic and magnetic energies were equal, with E = 1, He — 0.32 and 
Hm — 0.09. Phase diagrams are presented comparing the same b(k) for Runs 1A and 
1 in Figure 1: (a) k = (1,0,0) and (b) k = (0, 1, 1). In Figure 2, the time evolution of 
energy in modes with k 2 = 1, 2, 3 and 4, as well as the total energy for all modes, for 
real Run 1 is shown. 

Although the ensemble prediction is that Fourier modes such as b(k) are zero-mean 
random variables, this does not match the time behavior shown in Figure 1, particularly 
for k 2 = 1. Non-ergodicity is clearly evident and is due to the aforementioned broken 
symmetry. Also note that Figure 1 shows that j(k) = z'k x b(k) ~ b(k), which appears 
to be a general result, indicating ‘depression of nonlinearity’ [Kraichnan 1988], i.e., 
tendency to a force-free state where the Forentz force j x b ~ 0. Furthermore, the k 2 = 1 
modes track each other fairly closely initially; in Figure 1 (a), the ideal trajectory arrives 
at its attractor at t « 70, while the real trajectory peaks at l ^ 52 (the phase path ends at 
t = 300). The u(k) for k 2 = 1 also show similar behavior. 

Next, consider Figure 2, where the transfer of energy between shells of different k 
is shown. Although the total energy is monotonically decaying, the energy in k 2 — 4, 
3, 2 and 1 modes peaks at f ~ 2, 10, 20 and 50, respectively. This indicates an inverse 
cascade [Frisch 1975]; coupled with the attractor evident in Figure 1, we have a coherent 
structure or magnetic dynamo. 




CONCLUSION 


Although the simulations discussed were on a relatively small grid of 64 3 points, this 
allowed for runs from t = 0 to 300 (there is always this trade-off). In the ideal Run 1A, 
broken symmetry and non-ergodicity (i.e., ‘broken ergodicity’ [Palmer 1982]) leads to 
coherent structure similar to that found in recent 32 3 runs [Shebalin 2006] that evolved 
over t = 0 to 2000,. In the real Run 1, we saw a phase point trajectory close to that of 
ideal Run 1 A, up until t ~ 50. We expect that lower dissipation runs on larger grids will 
show even longer lived and more energetic magnetic dynamos. 

I would like to thank Dr. David Garrison of the University of Houston-Clear Lake for 
use of the Physics Program computing facility. 
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